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^ ! Abstract 

We present details for a method to compute the broken phase sphaleron rate (rate of 

^ ' hot baryon number violation below the electroweak phase transition) nonperturbatively, 

' using a combination of multicanonical and real time lattice techniques. The calculation 

•/^ ■ includes the "dynamical prefactor," which accounts for prompt recrossings of the sphaleron 

^ ! barrier. The prefactor depends on the hard thermal loops, getting smaller with increasing 

>• \ Debye mass; but for realistic Debye masses the effect is not large. The baryon number 

^ \ erasure rate in the broken phase is slower than a perturbative estimate by about exp(— 3.6). 

Csl ' Assuming the electroweak phase transition has enough latent heat to reheat the universe to 

^ ■ the equilibrium temperature, baryon number is preserved after the phase transition if the 

(X5 ' ratio of ( "dimensionally reduced" thermal) scalar to gauge couplings X/g"^ is less than .037. 



Ph! 1 Introduction 

^ ■ 20 years ago, t'Hooft showed that baryon number is not a good quantum number in the stan- 
dard model |]T|. The reason involves the nontrivial vacuum structure of the SU(2) (weak) 
gauge group of the standard model. In any gauge theory, the vacuum is not unique; any 
^ \ gauge transformation of A = has zero energy and is an acceptable vacuum. But SU(2) 
o3 ' (and any simple gauge group) has the property that the space of 3-D gauge transforma- 

tions is topologically nontrivial. A gauge transformation has an integer tts winding number 
associated with it. Since the winding number must be an integer, the space of smooth 
gauges, and also the space of vacua, is disconnected. The different connected components 
are characterized by their values of Chern-Simons number. 



which is an integer for a vacuum configuration, though not necessarily for an excited state. 
Classically, for the gauge fields to change from one topological vacuum at time t = ti to 
another at time t = tf, they must pass through excited states in the intervening time; to be 
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specific, 
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dNcs/dt = / d'xF^^F^'' , (2) 

which is clearly a gauge invariant quantity (though its integral, A^'cs, is not, because of the 
constant of integration). This means that it is possible to pass from a vacuum configuration 
to a gauge copy of that configuration via a path which cannot be smoothly deformed to 
remain always in vacuum. If we mod out the space of 3-D configurations by the gauge 
transformations, the space we get will then have noncontractible loops. 

This topological structure appears both in SU(2) (weak) and SU(3) (strong), where it is 
responsible for the physics of spontaneous chiral symmetry breaking. What t'Hooft noticed 
is that, because fermions couple to the weak SU(2) group of the standard model chirally, the 
anomaly relates Nqs to baryon number. If the gauge fields pass through some nonvacuum 
intermediate state from one topological vacuum to another (or around a noncontractible 
loop, if we think of configurations modulo gauge transformations), baryon number changes. 
Such changes are classically forbidden at zero temperature, so they only occur via quantum 
tunneling. Because the SU(2) gauge coupling is weak, and because the Higgs field breaks 
the symmetry, such processes are steeply exponentially suppressed, by ~ exp{—16iT^ / g'^) ~ 
10~^™. Hence such processes are of no terrestrial phenomenological interest. 

However, as a general rule, if a process only occurs in vacuum via quantum tunneling, 
then above some temperature it occurs much faster via thermal activation. (Chemistry and 
condensed matter physics are full of examples; annealing of crystal defects, for instance.) 
The same is true for baryon number changing processes in the standard model, although 
the "annealing temperature" is O(lOOGeV). In fact, there is a phase transition at Tc ~ 
lOOGeV in which the Higgs field loses its condensate, and above this transition baryon 
number violation is efficient |@, |^. It is quite possible that the baryon number of the universe 
originated in a cosmological electroweak phase transition, and in recent years this belief has 
driven the study of the electroweak phase transition. 

It takes two things for a cosmological electroweak phase transition to generate the baryon 
asymmetry of the universe. First, there needs to be enough CP violation to generate at least 
the observed abundance of baryons during the transition. Second, there cannot be too 
much "annealing" of the baryon number after the phase transition, that is, baryon number 
violation must be inefficient enough after the phase transition that a good fraction of the 
baryons survive to the present day. This is a condition on the phase transition's strength. 

The minimal standard model fails both conditions badly 0, ^ . However, well motivated 
extensions, like the MSSM with a light right scalar top, appear to be viable. Recent studies 
of baryon number production during the phase transition appear to show that enough CP 
violation can hide in places with few low energy consequences to generate the observed 
abundance of baryon number, and maybe a little more (for recent work see, for instance, 
1^, 1^, §). And the phase transition can be stronger. If the lightest scalar top is not very 
light, then perturbation theory can reliably relate the phase transition in the MSSM to the 
phase transition in the same effective theory used to study the standard model which 
has been well analyzed numerically 0, |1T|. If the lightest scalar top is lighter still, the 



phase transition may be stronger and more exotic |T^. This system can also be studied by 
nonperturbative lattice techniques [0. 

At present , the weakest link in our knowledge of baryon erasure after the phase transition 
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is the relation between the strength of the phase transition, now known nonperturbatively, 
and the efficiency of baryon number violation after it is over, for which we have only a one 
loop calculation 0, 0, We know that the perturbation expansion at high temperature 
near the electroweak phase transition cannot be viewed as an expansion in but at best 
as an expansion in a ratio of couplings. We also know that the two loop corrections in 
the perturbative expansion for the strength of the phase transition aren't very small in the 
"interesting" range of couplings where the baryon number violation after the transition is 
close to the efficiency limit. So it would be nice to actually know how good the one loop 
calculation is, or to replace it with a fully nonperturbative investigation. 

Very recently we have proposed a nonperturbative method to determine the rate of baryon 
number violation in the broken electroweak phase |T6[. This paper will ffil in all the details 
left out in that paper. Also, the calculation there was incomplete; it did not include a 
measurement of the "dynamical prefactor," discussed below. This paper will complete this 
aspect of the calculation. It will also discuss the importance to the broken phase sphaleron 
rate of hard thermal loops, which can modify the dynamical prefactor. 

For the impatient reader, we will present the basic ideas and the results right now. To 
find the sphaleron rate nonperturbatively, we ffist define nonperturbatively a surface called 
the separatrix, which sits half way between distinct topological vacua. Sphaleron transitions 
which permanently change iVcs must cross this surface. To find the Nqq, diffusion rate, we 
ffist compute the probability in the canonical ensemble to be in a narrow band about the 
separatrix; then we compute the mean inverse time to cross the band. The product is the 
probability flux across the separatrix. Then we compute a "dynamical prefactor," which 
tells what fraction of crossings lead to permanent resettling about a different vacuum. All 
three quantities can be computed nonperturbatively on the lattice, using a combination of 
Monte-Carlo and real time techniques. Including strong hard thermal loop effects modifies 
the dynamics in a way which lowers the dynamical prefactor, but for realistic parameter 
values the effect is minor. The iVcs diffusion constant is presented, and compared to an 
analytic estimate based on the two loop effective potential, in Table |I]. 

The paper is structured as follows. In Section |^, we outline the general idea of the 
calculation. Section |^ defines Chern-Simons number on the lattice, and the order parameter 
we will use, which is very closely related. It also discusses application of the definition to 
the symmetric phase case. Section | tells how we go about things numerically. Section ^ 
presents numerical results and compares them to a "semi-two loop" analytic estimate, and 
to the erasure bound. The last section concludes. For readers who are allergic to details of 
numerical studies, we recommend reading Section carefully, and perhaps the ffist subsection 
of Section ^, and then skipping to Section |. 

2 Broken phase measurement: general idea 

We want a technique for determining the A^^cs diffusion constant in the broken electroweak 
phase, where the rate is extremely small. The technique will be geared around the smallness 
of the rate and the fact that the system in finite volume will spend almost all of its time 
in a "neighborhood" of a topological vacuum, in a sense to be made precise below. These 
assumptions can be checked a posteriori, and do not constitute a real limit to the technique 
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1.38 ±0.02 


1.60 ±0.01 


1.82 ±0.03 


nonperturbative — ln(r(iT~^) (excl. prefactor) 


24.7 ±0.4 


28.3 ±0.4 


31.2 ±0.6 


— \n.{VdT~'^) (incl. prefactor) 


25.9 ±0.5 


29.5 ±0.5 


32.4 ±0.7 



Table 1: Perturbation theory versus nonperturbative F^. Appearances of T~'^ are really 
(2.5(?^Tc)~'^. X = \/g^ is the ratio of the Higgs self-coupling to the gauge coupling, and 
y = m'\j{T) / g'^T'^ is the dimensionless Higgs mass squared. The error bars for the nonper- 
turbative 00 are dominated by statistical errors in the determined value of Tc; errors in the 
nonperturbative value of F are statistical errors from the Monte-Carlo. The last row in- 
cludes the nonunity dynamical prefactor in the rate and should be taken as our most reliable 
estimate. 



in the broken phase. They will fail in the symmetric phase or when the phase transition is 
very weak, but in that case we can apply real time techniques 0, |18|, 0, |20|, ^ p3| , 
which can now produce quantitative results |^ (We should mention here that the 



symmetric phase case is not completely settled; it has recently been argued that there are 
logarithmic corrections to the parametric scaling behavior |]2^, which are however too small 



to be seen over the noise and other systematics still present in p6[| . 



2.1 thermodynamic approximations 

Before we start to describe our approach to the calculation we will specify the approximations 
to be made. 

We treat the thermodynamics of the standard model, or whatever extension is of interest, 

that is, as being well approximated 



in the dimensional reduction approximation |28 



by a three dimensional, bosonic path integral with parameters carefully matched to those of 
the full theory. This is an excellent approximation and we have no regrets in making it. For 
a study of corrections to this approximation in the present context, see [^|, which shows 



that the leading thermodynamic effects not included in the dimensional reduction procedure 
have a negligibly small effect on the sphaleron energy. 

Conveniently, dimensional reduction is equivalent to treating the theory's thermodynam- 
ics as equivalent to those of the classical bosonic theory [0, with certain mass counterterms. 
Similarly, we can treat the theory's dynamics in a classical approximation. This should be 
valid in the infrared 
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j], with one serious complication. That is, the structure of di- 
vergent radiative corrections to unequal time correlators is much more complicated than 
that for equal time correlators. For the equal time correlators, which are all that matter 
to thermodynamics, the divergent radiative corrections are mass squared corrections for the 
Higgs and Aq fields, which can be computed once and balanced by counterterms. For un- 
equal times the linearly divergent radiative corrections, the hard thermal loops, have a more 
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complicated structure. And they can significantly change the dynamical behavior on time 
scales longer than the inverse plasma frequency. The modifications can be important for the 



Chern-Simons number diffusion rate |33| , |34| , [35| , as has recently been verified numeri- 



cally for the symmetric phase dynamics In the current context they will modify the 

"dynamical prefactor", but they have little bearing on those parts of the calculation which 
are thermodynamical. 

We will also frequently make the approximation that the so called "heavy" degrees of 
freedom can be integrated out |^ , including both the time component of the gauge fields and 
any squarks present, so the theory reduces to an effective theory for the minimal standard 
model 0. This theory is specified by two parameters; x = X/g"^ the ratio of scalar to gauge 
self-couplings, and m^(T) the thermal Higgs mass squared Lagrangian parameter. m^(T) 
is a monotone increasing function of T, going from 0{g^T'^) at high temperatures to quite 
negative (symmetry breaking) in vacuum. The phase transition occurs near where it is zero, 
at a critical temperature which in the context of dimensional reduction becomes a critical 
m^iT). At tree level, X/g"^ = {mn/my^y /8 but the radiative corrections are important and 
the relation between X/g"^ and the ratio of physical zero temperature masses mn/mw is not 
simple, especially in extensions to the standard model with new light bosons. Also note that, 
for instance, if the MSSM right scalar top is too light, the reduction to an MSM like effective 
theory is not very reliable; we should use an effective theory which contains the light squark 
and the gluons. Dropping heavy modes is not a necessary step for using our technique, it 
is merely convenient to reduce the numerical demands, which would be a few times larger 
if we include the Aq fields and order 10 times larger if we include the squarks and QCD. 
We will discuss how we think light squarks would change our results in the conclusion. We 
discuss the matter of integrating out the Aq field in more detail in subsection |4.2| . It is not 
always appropriate to do so, and in particular we cannot when we are studying the influence 
of hard thermal loops. 

We do keep the U{1) subgroup, which is often left out in electroweak studies. Its role in 
setting the sphaleron rate is probably almost entirely due to its effect on the strength of the 
phase transition and not a direct modification of the sphaleron, see [^, but the numerical 
cost of including it is small enough that dropping it is pointless. We use tan^ 9w = -32, 
based on a 1 loop match between vacuum MS and 3-D thermal values using results in [RO . 



2.2 the separatrix 

The idea of the separatrix between vacuum states is essential to our technique. Before 
introducing it, let us review what we expect the space of gauge-Higgs configurations to 
look like in the broken phase. The space of three dimensional gauge-Higgs configurations 
is periodic, with a discrete set of vacua. To be more precise, we should consider the space 
of gauge-Higgs configurations modulo (all) gauge transformations. In this case all vacua 
coincide^, but the space is not simply connected. Since the index of the Dirac operator 

^If the global topology of space is multiply connected then Yang-Mills theory has a connected manifold of 
inequivalent vacua corresponding to different values for traces of certain noncontractible Wilson loops. The 
toroidal spaces we will consider are multiply connected; however, the (fundamental representation) Higgs 
condensate lifts the degeneracy of the would be gauge vacua. These complications will not be important for 
what we do. 
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Figure 1: Cartoon of periodic vacuum structure, separatrix, typical path which stays near 
a vacuum, and exceptional path which crosses the separatrix and leads to permanent Nqs 
change 



remembers when we go around a noncontractible loop in configuration space, the relevant 
space of physical configurations is the universal cover of the space of configurations modulo 
(all) gauge transformations. The universal cover has a discrete set of vacua labeled by 
the index of the Dirac operator^. If the line connecting two vacua in the universal cover 
projects to a winding 1 loop-or in more conventional language, if two vacua differ by 1 in 
Chern-Simons number-we will refer to them as neighboring vacua. 

We expect that, in the broken phase, almost all of the weight of the canonical ensemble 
lies in states which are in some sense close to one of the vacua. The Hamiltonian evolution 
of a generic state in the ensemble will wander around in the neighborhood of one vacuum 
for an exponentially long time before it happens to make an excursion far enough away that 
it crosses to being nearest another vacuum, which it may find instead. We expect that this 
is how A^cs diffusion will occur. 

To determine the rate of A^^cs diffusion, we draw a (codimension 1) surface separating 
one vacuum from its neighbor, halfway between the minima in some sense, so that all of the 
well populated area near one minimum falls on one side and all the well populated area near 
the other minimum lies on the other side. This surface is called the separatrix between the 
vacua. To cross from being near one vacuum to being near another vacuum, a Hamiltonian 

■^The universal cover is roughly the same as the space of configurations modulo small gauge transforma- 
tions. But we prefer to think in terms of the universal cover of the space of configurations modulo all gauge 
transformations, because this is an explicitly gauge invariant approach. 
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^ ^ "Good" separatrix 




Figure 2: Cartoon of how a poor choice of separatrix can lead to overcounting the flux, and 
a small dynamical pref actor. 



trajectory must pass through the separatrix dividing them. Wc will assume that, after such a 
crossing, the trajectory almost never promptly continues to and crosses the next separatrix, 
but instead either settles around the new minimum for long enough that ergodicity "erases 
its memory," or turns around and returns to the vacuum it started from. Then the flux of 
probability of the thermal ensemble through the separatrix is an upper bound on the diffusion 
rate for Nqs- It is an upper bound because of trajectories which cross the separatrix, turn 
around, and return to the original vacuum. These lead to flux of probability through the 
separatrix, but not to Ncs diffusion. To get the true diffusion rate, we need to find not 
only the flux through the separatrix, but the average number of crossings of the separatrix 
per long term change from being near one vacuum to being near another. We will call the 
reciprocal of this, the fraction of separatrix crossings which are associated with permanent 
Nqs change, the "dynamical prefactor". If we know both the flux of probability across 
separatrices, and the dynamical prefactor, then we know the diffusion constant for A^cs; the 
diffusion constant is 

^ ((iVcsW - iVcs(O))-) ^ , ,3, 

t— >oo t 

What we really want is the diffusion constant per unit volume, = 7d/V^. We should use a 
volume which is large enough to prevent finite size systematics but small enough that there 
are almost never two simultaneous A^"cs changing events in different places. 

A few points are in order here. First, the assumption that there are no prompt crossings 
of multiple separatrices is essential to the calculation. As we will see, it is also easy to test. 
Second, if we make a somewhat poor choice of the separatrix, so that there is some place 
where it bends nearer to one vacuum than the other, then much of the fiux in that place will 
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be of trajectories which double cross and return to their starting vacuum. The probabihty 
flux will be larger than with a better definition of the separatrix. However, the dynamical 
prefactor will be smaller, since these extra crossings do not lead to permanent Nqs change. 
The rate we determine is independent of exactly where we put the separatrix as long as the 
flux across it is exponentially small, and as long as we make a complete calculation, including 
the dynamical prefactor. We illustrate this point in Figure |^. In practice we should look 
for a good choice of separatrix, since a poor choice of separatrix may make it harder to get 
good statistics for F, as most of the (numerical) effort will go towards studying trajectories 
which double cross rather than ones which really change A'^cs permanently. 



2.3 calculation: perturbation theory 

The existing perturbative calculations of the broken phase diffusion constant for A^^cs are 
along the lines of the approach we just described. To allow a perturbative calculation, they 
make an additional assumption; that the separatrix is dominated by its saddle point. That 
is, they assume that gauge-Higgs configurations on the separatrix look like perturbatively 
small excitations about a background field which is the lowest energy point on the separatrix, 
which will be Klinkhamer and Manton's sphaleron |3^. The choice of a definition of the 



separatrix is then made perturbatively; a point is on the separatrix if the excitation in the 
unstable direction of the sphaleron is zero. The probability flux through the sphaleron can 
be computed in perturbation theory by comparing the free energy of all excitations of the 
sphaleron to the free energy of excitations about the naive vacuum, with the frequency of 
the unstable mode serving to convert a probability into a flux and the translational zero 
modes converting this into a flux per unit volume. 

The probability flux through the separatrix has been computed in the above approxima- 



tion at the one loop level [|14[ |15| . However, extending the calculation beyond one loop raises 
severe technical problems. The sphaleron is not a spatially homogeneous background field, 
so the perturbative calculation must be done in real space with a numerically determined 
spectrum of fluctuations. The one loop calculation requires finding this spectrum, but the 
two loop calculation involves overlap integrals to compute the energy of their mutual in- 
teractions. There are also conceptual problems, because one of the fluctuation directions is 
unstable. At one loop it is excised from the sum over fluctuations; the other fluctuations set 
the probability to be near the separatrix and it turns that probability into a flux through 
the separatrix. But it is not clear how to separate it from the other modes at two loops. 
These problems obstruct a systematic improvement of the perturbative treatment. 

There is also the problem of how to determine the dynamical prefactor perturbatively. 



Khlebnikov and Shaposhnikov argued that it should equal 1 |^ , but Arnold and McLerran 
made an estimate based on Landau damping which suggests that it is quite a bit less than 1 
0]. That argument has been more carefully developed by Arnold, Son, and Yaffe, who claim 
that, so long as the Higgs condensate gives a W mass which is parametrically mw ~ g'^T, 
the prefactor should be parametrically 0(q;^) [Q. Their argument has recently been tested 



in the symmetric phase [^. No one has used their picture to get a quantitative prediction of 
the dynamical prefactor within the context of the perturbative calculation of the sphaleron 
rate, although this should be possible in principle. 

These limitations of the perturbative approach, together with the generally spotty per- 
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formance of perturbation theory for electroweak phenomena at temperatures near the elec- 
troweak phase transition, motivate a fully nonperturbative attack on the problem. 



2.4 topology and the lattice 

The nonperturbative technique best suited to studying the sphaleron rate is the lattice. 
There is an apparent complication, though, which is that topology cannot in general be 
well defined on the lattice; the global structure of the space of three dimensional lattice 
gauge-Higgs configurations is different than in the continuum, and in particular it is simply 
connected, so there are not topologically distinct vacua. However, if we make the lattice 
spacing suitably small, then the thermal ensemble is completely dominated by configurations 
in which every elementary plaquette is close to the identity. This means, roughly, that 
fields are perturbatively small at the lattice scale. This subspace of the space of lattice 
configurations does have the same topological structure as the continuum theory. Hence, if 
we excise a subspace of lattice configurations which carries an exponentially small weight, 
then we can talk about topology on the lattice. 

The physical meaning of this is as follows. The space of continuum configurations permits 
sphaleron-like objects of arbitrarily small spatial extent. They also exist on the lattice, down 
to where their size is comparable to the lattice spacing; but at this point it becomes unclear 
how to define a smoothly interpolating continuum field, and the topological meaning is lost. 
Such a configuration in the analogous 4 dimensional context is referred to as an "exceptional 
configuration." However, the energy of spatially small, sphaleron-like field configurations 
rises linearly with inverse radius; so the Boltzmann suppression of lattice scale sphalerons 
is enormous and they essentially never occur. As long as the "genuine" sphalerons we 
study are comfortably larger than the lattice spacing, we have no problem. And when the 
genuine sphalerons we want to study are not comfortably larger than the lattice spacing, 
then obviously the lattice spacing is too coarse, and we should use a finer lattice. 

The situation here is much better than it is in the 4-D case considered in QCD. It is also 
true in 4 dimensions that if the gauge fields are smooth enough, then topology is well defined 
However, in practice the fields may not generally be smooth enough. Instantons 



just larger than the lattice spacing typically do exist, because the instanton action is classi- 
cally scale independent, and so only varies logarithmically with instanton size. Making the 
lattice finer does not eliminate lattice spacing sized instantons (exceptional configurations) 
very quickly; in fact, because their density goes as exp{—l/g^) and g"^ varies logarithmically 
with lattice spacing, their density declines as an algebraic power of a. In the 3-D context, 
though, the energy of a sphaleron of radius r goes as l/{g'^r) and the density of exceptional 
configurations varies with lattice spacing as exp(— (coefficient)/((7^aT)), which falls off ex- 
tremely rapidly as a is made small. This difference between the 3-D and 4-D cases is because 
the 3-D theory is super-renormalizable; the coupling constant is g'^T, which is dimensionful. 
Since l/g"^ appears in the exponent for the rate of any nonperturbative phenomenon, and 
since it must be accompanied by a T, on dimensional grounds a nonperturbative phenomenon 



*in fact the 3-D case is a subset of the 4-D one, since the topology we are talking about is the second Chern 
class of a closed loop in 3-D configuration space, which is equivalent to a periodic 4-D lattice configuration 
with the spacing in the fourth dimension driven to zero. 
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which occurs at the lattice spacing scale must proceed at a rate which goes as 1/a, leading 
to the exponentially fast rolloff in the density of exceptional configurations as a ^ 0. 



2.5 a "nice" nonperturbative choice for the separatrix 

Now that we have established that topological questions can have a meaning on the lattice, 
we ask how to choose a separatrix which is defined nonperturbatively and can be implemented 
on the lattice. 

A definition we would prefer for the separatrix is the "gradient flow" definition, suggested 



in [^, ^ . It is defined in terms of gradient flow under the Hamiltonian. The Hamiltonian 
is a smooth function over the space of configurations, with degenerate global minima at the 
vacua. We believe that these are the only local minima of the Hamiltonian in 3-D Yang- 
Mills Higgs theory, although we do not know a proof. Hence, following the direction of 
steepest descent of the energy (gradient flow) will lead, off a set of measure zero, to a vacuum 
configuration. A rigorous definition of "the vacuum closest to a configuration" is the vacuum 
arrived at by such gradient flow, which is also easily implemented on the lattice. A very 
sensible definition of the separatrix is then the boundary between the gradient flow basins of 
attraction of two neighboring vacua (vacua with A^^cs differing by 1). Equivalently, it is the 
basin of attraction of the saddle point which sits between the two vacua, i.e. the sphaleron. If 
we mod out by all gauge transformations the two vacua are equivalent, but the separatrix can 
still be defined as the surface where two infinitesimally separated configurations on opposite 
sides will have macroscopically different gradient flow paths which, when spliced together, 
form a noncontractible loop. 

Alternatively we could define the separatrix just in terms of the Yang- Mills field (connec- 
tion) and the Yang- Mills term in the Hamiltonian. This choice has the added benefit that, 
as the configuration gradient flows under the Yang-Mills Hamiltonian in 3 dimensions, it 
becomes exceedingly smooth (meaning that all gauge invariant local measurables are slowly 
varying and the energy density is very small). Also, leaving out the Higgs fields evades the 
complication that the Higgs mass squared is renormalized by ultraviolet thermal excitations, 
which change during the gradient flow. The Yang-Mills gradient flow separatrix may not 
coincide with the Yang-Mills Higgs separatrix (which is not defined until we decide how to 
deal with the renormalization of the Higgs mass term in the Lagrangian). But if we perform 
a complete calculation, including the dynamical prefactor, then the exact choice of separatrix 
should not matter, as long as crossings are exponentially rare. Of course, a poor choice will 
make the calculation inefficient, since most crossings will not be associated with topology 
change. But we do not expect the Yang-Mills gradient flow separatrix to be a poor choice, 
and we will be able to check this belief when we study the dynamical prefactor. 



2.6 approach to computation 

Now, we will outline how to compute the flux and the dynamical prefactor. To do so, we 
need more than just a definition of a separatrix. We need an order parameter which takes 
a special value, say A^ = 1/2, on the separatrix, and is smaller on one side and larger on the 
other (say, going from A^ = at one vacuum to A^ = 1 at the other). Assume that we have 
such an A^. 
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First, we should measure the probabihty over the canonical ensemble that is within 
some small tolerance e/2 of (integer+)l/2. This gives the probability to be very near a sep- 
aratrix. Since the probabihty to be close to the separatrix is (expected to be) exponentially 
small, we will need to use a multicanonical reweighting technique to sample here accurately. 
The basic idea is presented early in Section |[ 

Then we need to know the mean inverse time for crossing this narrow region, to turn 
the probability into a flux. This is 1/e times {\dN/dt\), the mean of the absolute value of 
the time derivative of N , where the averaging is over the ensemble restricted to the narrow 
band about the separatrix. We determine this by taking a canonical sample of states with 
|A^ — 1/2| < e/2, drawing momenta for each from the thermal ensembleQ, and performing 
a short segment of Hamiltonian evolution, comparing with its value a short Hamiltonian 
time later. Thus, we can determine the flux. 

To determine the dynamical prefactor, we take a canonically weighted sample of configu- 
rations at the separatrix and choose momenta for each, just as we do to find {\dN/dt\). Then 
we perform the Hamiltonian evolution both forward and backward in time, until the forward 
and backward histories both settle into the neighborhood of a vacuum. We count how many 
times the Hamiltonian trajectory crosses the separatrix before it settles semipermanently 
about a minimum. The prefactor is 



1 x<^ final vacuum 7^ initial vacuum 

number of crossings 1 final vacuum = initial vacuum 



sample 



(4) 



This is not the same as adding up the number of times the final vacuum differed from the 
initial one and dividing by the number of crossings of the separatrix. The reason is that the 
latter overcounts trajectories with many crossings, since the ensemble samples them more 
often than trajectories with fewer crossings. 

This is our recipe; given an order parameter N , we can determine both the flux of 
probability through the separatrix, and the dynamical prefactor, and hence F^. 

It remains to choose an order parameter. We want one such that the N = 1/2 separatrix 
will be either the same as or quite close to the "good" choice of the gradient flow separatrix. 
Naively, we expect a good N to be the Chern-Simons number, Nq^- Actually, this is not 
quite right, as we discuss at some length in the next section. 



3 Defining lattice A'^cs 

We now have an idea for a nonperturbative definition of the separatrix, but we need an 
order parameter which varies from to 1 as it ranges between vacua and equals 1/2 on a 
surface close to the Yang-Mills gradient separatrix. The literature generally considers the 
Chern-Simons number Nqs to be the best choice for this so we will discuss how to 



define Nqs on the lattice; then we will define a different but closely related order parameter 
A^, which gives a separatrix much closer to the gradient flow separatrix than A'^cs does, and 
which we will actually use in the calculation of the broken phase diffusion rate. 

^It is essential here that the phase space is the tangent bundle of the configuration space, and that the 
thermal probability distribution is a product of a configuration space function and a (Gaussian) function 
over momentum space. 
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3.1 the definition 

Chern-Simons number should be defined as some real valued function over the space of three 
dimensional gauge connections. Let us review some properties which A^^cs would have in the 
continuum, and see how closely we can preserve them on the lattice. 

In the continuum, Nqs, defined on the space of gauge field configurations modulo all 
gauge transformations, has the following properties: 

1. Ncs should be a continuous, multiply valued function, with the values separated by 1. 

2. The multiple values correspond to the projections of different points on the universal 
cover, where A'^cs is single valued. A noncontractible loop in the configuration space lifts 
to a line segment in the cover space, and the difference in A'cs between the endpoints is 
the winding number of the loop. (A more conventional way of saying this is that A^cs 
differs between two gauge copies by the winding number of the gauge transformation 
between them. But we prefer the above, gauge invariant statement.) 

3. A vacuum configuration has A'cs modulo 1 equal zero. 

4. Consider a path in configuration space parameterized by r. The difference in A'cs 
between the beginning and end of the path should be 

AiVcs = ^l^^l d'xe^.kFt.iD^AkT • (5) 
This is the same as saying that A"cs is the indefinite integral of 

d'^F'^J;^ = d'.EtB; . (6) 



How many of these properties can we preserve on the lattice? Not all of them; as we 
pointed out in [20|, there is no local operator definition of -Ef-Bf which is a total derivative. 
We also argue there that there are severe difficulties satisfying 2.; but this is because we 
demanded a singly valued, non gauge invariant definition of A^cs- In fact we can present a 
lattice definition of A^cs which satisfies everything but 4. and the continuity requirement, 
provided we restrict to gauge fields which are smooth enough, in the sense of the last section, 
so that 2. makes sense. 

Before constructing a lattice definition of A'cs, we remind the reader how the lattice fields 
are defined. On a lattice, scalar fields are only defined at a discrete set of points, the lattice 
sites. The gauge field should be a connection, that is, it should be a rule which tells how 
to parallel transport fields along paths. We allow a path on the lattice to consist of a series 
of straight lines between nearest neighbor lattice sites. The connection is then defined by 
associating a group element U gSU(2) with each of these elementary straight lines (referred 
to as the links of the lattice). A small closed path, or the product of the U around the path, 
is called a plaquette; the 1x1 square is the elementary plaquette. When we do not specify 
the shape of a plaquette we mean an elementary plaquette. The product of the U around 
a plaquette (written [/□, often just referred to as "a plaquette") will not in general be the 
identity; its failure, a curvature in the connection, is a field strength. 
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Figure 3: Illustration of the pieces which make up a lattice gauge theory. 



It is important that the field strength is not associated with a site of the lattice, but 
with a plaquette, which sits in between sites. Then, to make a lattice implementation 
of / e^'^°'^ we will have to do some averaging. The argument of the integral is a 
pseudoscalar and should perhaps be defined at the lattice sites; but the field strength is 
associated with a plaquette. To preserve cubic symmetry, F^^, at a site will have to be the 
average over the four plaquettes in the /x, v plane which touch the site.0 Because of this 
averaging process over things which do not live quite at a lattice site, the resulting lattice 
definition of e^^'^^ F^j^^Fap will not be a total derivative; and we cannot fix this problem by 
going to fancier definitions involving weighted averages of plaquettes of various shapes [pO |. 
The problem is that the continuum proof that FF is a total derivative relied on continuity 
of the fields, and this continuity is lost on the lattice. This makes it impossible to satisfy 4. 
Note, however, that the lattice definition of F^^F^^ is gauge invariant. 

However, when the gauge fields are weak at the lattice scale (meaning that the elementary 
plaquettes are close to the identity) and slowly varying (meaning that the departure from 
the identity is nearly the same between a plaquette and the parallel transport of a nearby 
plaquette in the same plane), then the lattice definition of E^B"^ is approximately a total 
derivative. Further, we can pick a definition of iVcs so that 4. is satisfied at least in one 
particular special direction. We choose it to be true in the cooling direction, that is, the 
direction of energy gradient fiow. Here the energy of the lattice field is given by the standard 
Kogut-Susskind Hamiltonian which is defined up to a multiplicative factor as 



plaquettes Ud 



(7) 



The gradient of i^KS is to be understood in terms of the metric of the configuration space, 

^We are thinking of a path through configuration space as a 4 dimensional lattice, where the fourth 
dimension is the path parameter, taken as a discrete rather than a continuous variable. Successive 3 di- 
mensional slices are configurations along the path. In a numerical setting, a path through the configuration 
space will always look like this, though making the spacing in the 4'th direction arbitrarily small restores 
the continuity of the path. 
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Figure 4: Original lattice and the coarsened version, which is just the solid lines and filled 
vertices. The link matrix between two vertices of the coarsened lattice is the product of the 
matrices on the two original links which make it up. 



which is the product of the Haar measure over each link matrix U. 

We will call a path which follows the steepest descent of i^KS a cooling path, and we 
parameterize it with a cooling time r, defined as dr = (i(path length) /[rfi^Ks/c^lpath length)]. 
The gauge fields evolve along the path according to |25| 



^ = -D^UD'^Hks , (8) 

where is the left acting covariant derivative, D°'U = ir'^U, and D^'Hyis is i^KS with U 
replaced with D°'U. This is the gauge invariant lattice implementation of the continuum 
evolution 

dAtjx) ^ dH 

dr dA^{x) ■ ^ ' 

T has dimensions of length squared. 

We are defining A^cs so that condition 4. is true along the cooling path. Since cooling 
eventually leads to the vacuum off a set of measure zero, and since A^cs of the vacuum is by 
definition an integer, we get 

Ncs = integer - H dr^ [ d^xE^Bf , (10) 
where we mean the lattice definition of EfBf which we have described, and which is written 



down in [|T9[. This definition of A^'cs is nice because the cooling path leads most quickly 
towards configurations with weak, slowly varying gauge fields, so the definition is minimally 
contaminated by the problems with the lattice definition of E^B^ which we have discussed. 

Performing the integral in Eq. (|To|) takes an enormous numerical effort. Cooling the fields 
(following the cooling path) stably requires using a r step size of less than (a the lattice 
spacing), or else the most ultraviolet excitations become unstable.[] But sufficient cooling 

^We improve on this marginally by using alternately larger and smaller step sizes, analogous to what we 
did in pol when quenching the E fields to enforce Gauss' Law. 
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Figure 5: {g'^/Sir'^) J EfBfd^xdr for a series of points on a Hamiltonian trajectory. It is 
clear where to adjust the integer part of Eq. ([T0|) to keep Nqs (approximately) continuous. 

may demand going to r of order lOOOa^. However, the cooling is by far the most efficient at 
removing ultraviolet excitations, and already by r > the fields are slowly varying. After 
this much cooling we lose almost no information if we drop some UV degrees of freedom 
by setting up a coarsened lattice. Define an even site as a site where all three coordinates 
are even numbers. In a scalar field theory, we would coarsen by dropping out all the lattice 
sites which are not even, leaving a lattice half as many points across in each direction. In 
a gauge theory, we also need to define the connections between the sites of the coarsened 
lattice; we define the connection between two neighboring even sites as the product of the 
two connections along the straight line between them. We illustrate the idea in Figure ^ 
The remainder of the cooling then proceeds 2^ times faster, 2^ because the lattice is smaller 
and 2^ because we can use a r step size which is larger in physical units. 

Of course we must check that this procedure produces the same answer as we get by not 
coarsening. But if we use an O(a^) improved definition oi E ■ B, for instance the one we 
developed in |^], then this is not a problem at all. The value of Ncs gets rescaled by an 
amount typically less than 1% and has an amount of noise added to it which is typically even 
smaller. For lattices more than about 28 sites across, we can even safely perform a second 
stage of coarsening after performing several coolings on the once coarsened configuration. 



3.2 application to the symmetric phase 

As an application we discuss how to use this definition of A'^cs to study A'cs diffusion in the 
symmetric phase, or pure Yang-Mills theory. In this case we want to track A^cs along the 
projection into configuration space of a Hamiltonian trajectory in phase space. That is, gen- 
erating a Hamiltonian trajectory will give us a closely spaced series of points in configuration 
space, and we must associate a single valued function A'cs with this series. Our definition 
says that we should measure A"cs at each point by performing the integral in Eq. ([T0|), 
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Figure 6: (Time, cooling time) plane, with curves used to track A'^cs for a Hamiltonian 
trajectory. Every few Hamiltonian updates, we construct a cooled copy of the configuration, 
giving a parallel cooled path. Every few points on that path, we measure iVcs directly, and 
we fill in between by integrating E ■ B along the parallel cooled path. As long as the integrals 
along paths 1 and 2 always agree modulo an integer, within a small tolerance, the technique 
is topological. 



that is we must cool the configuration at each point along the Hamiltonian trajectory to the 
vacuum. We also have to choose the integer part of Eq. (^) somehow. At certain points, 
the value of the integral will abruptly jump by almost an integer. This happens whenever 
the Hamiltonian trajectory crosses a gradient fiow separatrix. We should choose the integer 
part of A'cs to minimize the magnitude of the jump in A'cs at the separatrix. We illustrate 
with an example of real data from a simulation of Yang-Mills theory in Fig. ^. Note that 
the discontinuities in Nqq, do not necessarily occur when A'cs is near ±1/2; we will discuss 
this more in the next subsection. 

In practice, even with lattice coarsening, the numerical costs of cooling every configuration 
are unbearable, but we can do better. We show how in Fig. Every few steps, we cool 
a little, to a depth of r ~ a^, and we thereby construct a cooled image of the Hamiltonian 



path, a technique explored by Ambj0rn and Krasnitz p5|. We measure A'cs using Eq. (|10[) 
at occasional points along this path, interpolating A"cs in between by integrating E°'B°' 
along the cooled image of the Hamiltonian path. The cooling has eliminated most of the 
UV excitations, so EfBf along the cooled path is close to behaving as a total derivative, 
especially if we use an O(a^) improved definition of E^Bf. The interpolated value of Nqq, is 
therefore almost what we would get by using the direct definition at each time. We only need 
it to be close enough to determine the integer part of Eq. ( ]TI1| ) unambiguously, which we can 
if the value we get by integrating on path 1 in Fig. ^ and the value we get by integrating 
on path 2 in that figure differ by an integer plus a small remainder. We can think of this 
remainder as a calibration of the integration along the cooled path, so the approach is a 
"calibrated cooling method" , with the occasional cooling paths to the vacuum recalibrating 
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Figure 7: Left: A^'cs measured by the slave field (upper curve) and "calibrated cooling 
method" (lower curve). Right: the difference between the curves (note scale). The difference 
is small and spectrally white, so the methods are in good agreement. 



the method of Ambj0rn and Krasnitz ||2^ to make it topological. 

We compared this approach of measuring A^^cs to the "slave field" topological method 
2^ , by evolving Yang-Mills theory on a 24^ grid at a = l/{2g'^T) (/5l = 8) for a total time of 



6000a, tracking Ncs by each technique. For the technique we just described, we constructed 
a cooled image path with one point every a/5 time. The cooling depth to this path was 
5a^/8, and we calibrated by cooling to the vacuum every 2a time. With this lattice spacing 
and this frequency of calibrating, the largest remainder we observed was 0.2 and the typical 
absolute value was less than 0.05. We present the results for Nqs in Fig. ^. We have offset 
A'^cs measured by the slave field method by 5 to keep the curves from lying on top of each 
other. The agreement is outstanding, and the difference in the determined values of Nqs is 
white on long time scales. 

To explain why the two methods have a white noise difference, we review briefly how the 
slave fleld method works. It tries to keep track of the "integer" part of Eq. (|10D by assuming 
that the cooling path will end in a vacuum which has winding number zero in Coulomb gauge, 
and then adding up the number of large gauge transformations required to keep the system 
in Coulomb gauge during the Hamiltonian trajectory. However, this ignores the contribution 
to Eq. (p!OD from the integral, that is, the difference between Nqs of the conflguration and 
of the vacuum it cools to. Also, the algorithm used to flnd Coulomb gauge sometimes 
gets trapped temporarily in a Gribov copy with a different winding number. But neither 
difference between the methods will grow without limit in time, so the difference between the 
two measurement methods for A'cs is white on long time scales. Hence, the derived diffusion 
constant will be the same within errors (caused by the white noise difference). Indeed, when 
we used the technique of to extract from each trajectory, the two methods of tracking 
A"cs gave the same answer within error (7^ = .0515 ± .0077a"^ for the new method, versus 
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7d = .0516 ± .0082a~^ for the slave method). 

As an aside, we mention that Hetrick and de Forcrand have used the method of coohng to 
resolve the Gribov problem, by defining Coulomb gauge (or, in their 4 dimensional context. 
Landau gauge) as the gauge in which the vacuum configuration arrived at via cooling has 
the minimum value of / 
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3.3 modifications for measuring in the broken phase 

The definition of Nq^ which we just presented more accurately reproduces the continuum 
meaning of A'^cs than any other we know. But we don't actually want all of the attributes 
of the continuum meaning of A^cs if "we want to measure Td in the broken phase by the 
separatrix method. The reason is that A"cs is only directly a measure of topology for vacuum 
configurations. There are contributions to A^cs in excited states which are uncorrelated to 
topology; for instance, A'cs does not vanish in abelian gauge theory, even though that theory 
has no tts topological sectors. In the continuum abelian theory, A"cs is given by 
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and the mean square value of A'cs is 
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Using Wick's theorem and the momentum representation of the propagator, in a general 
covariant gauge, this becomes 
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(13) 



so A"cs will be Gaussian distributed with a linearly divergent variance. On the lattice, the 
UV divergence will be cut off by the lattice scale; the coefficient was found by Amjorn and 



Krasnitz 19 and is 



(Ar2g) = (1.44xlO-V^TVa, 



(14) 



The divergence occurs because, while the energy cost of storing A^cs in a UV mode grows 
linearly with p, the number of available states in which to store A^cs grows faster; entropy 
wins over energy. The same thing happens in SU(2) theory, because A"cs also contains the 
eijkFijAk term. The coefficient of the divergence in SU(2) is larger by 3, the dimension of 
the group. Since Yang-Mills theory in 3-D is super-renormalizable, the UV decouples from 
the IR, where the genuine topology changing physics occurs, so this UV divergent, Gaussian 
contribution will appear as an additive correction to the IR contribution to A"cs. That is, to 
reasonable accuracy we can think of A'cs, defined in Eq. (0), as A'cs = A'^g + Nq^ , where 
topological information is in Nq^, and A'^g' is independent of A^^g and Gaussian distributed. 
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Given a single IR field configuration with some particular value of Nq^, different realiza- 
tions of the UV excitations on top of the IR fields will then give a distribution of values of 
Ncs', so the probability that Nqs will have a particular value x, Ptsscs{x), will be 

^Ncs(x) = / / 6{x- NEI - N^I) , (15) 

jIR configs J\JY excit. 

or, defining P^si^) and -Pncs(^) ^'^ ^e the probability distributions for the IR and UV 
components of A'^cs, 

P^csix) = jdyj dzP^c,iy)P^^,iz)6ix -y-z) = J dyP^%iy)P^^,ix - y) . (16) 

The probability distribution for A^^cs will be a convolution of the interesting IR distribution 
and a Gaussian noise distribution. Convolving any periodic function with a Gaussian always 
degrades contrasts in that periodic function, enhancing the probability to be near the A^^cs = 
1/2 separatrix relative to the case of no UV noise. The distortion will be small only if 



dx'^ 



(17) 

=1/2 



Later, we will use a 40^ lattice with a = 2/{5g'^T). Multiplying Eq. (0) by the group 
factor of 3 and plugging in numbers, ((A'^g')^) = 0.44 for such a lattice, which is too big. So 
defining A'cs by Eq. (p!0D will not do. 

Another way to state the above is that the separatrix one gets from the condition A'cs = 
1/2 is sensitive to ultraviolet excitations, which makes it "all wiggly"; it will have lots of 
"fingers" which stick out towards one or the other topological vacuum, and the problems we 
discussed in Section 0, of there being many crossings of the separatrix which don't have to 
do with permanent A"cs change, will be severe. 

The problem is that A'cs = 1/2 is not particularly similar to the "good" gradient flow 
definition of the separatrix. We want an order parameter which is close to A^ = 1/2 on the 
gradient flow separatrix. The N = 1/2 separatrix does not need to correspond exactly with 
the gradient flow deflnition; is sufficient if the distribution of values of A^ on the gradient 
flow separatrix is narrow, preferably narrower than d'^P^{x)/dx'^ . A slight change to Eq. 
(p!0|) will do the trick; deflne 

A^ = integer + ^ T dr [ d^xE^B^ , (18) 



d7r2 Jr, J 

where we mean the lattice implementation of the integrals and E, B as before. In other 
words, we "pre-cool" the conflguration for cooling time tq and then measure A"cs- The 
pre-cooling is intended to remove UV excitations without affecting the underlying IR flelds 
much. 

Now consider A^ in the abelian theory again. The theory is linear, so it is easy to analyze 
how cooling affects it. A particular transverse mode A{k) evolves according to 

dA{k) dH _ 



dT dA{k) 



-eA{k) A{k,T) = e-''^A{k,0). (19) 
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Figure 8: Ai'cs (wildly oscillating curve) and (curve which stays near 0) during a broken 
phase Hamiltonian trajectory, in a 40^ box with a = 2/{5g'^T) (and using tq = 3.2/((7^T)^). 
Nqs has a lot of UV noise, which is absent in N. 



The propagator in Landau gauge becomes 

{Ai{k,T)Aj{l,T)) = 



-6{k + l) 



and the variance of N is 



(20) 



(21) 



647r4 7 (27r)3 (p2)2 • 

So pre-cooling removes the UV noise from the definition of A^^cs- 

For comparison. Figure ^ shows N and A^^cs for the same broken phase Hamiltonian 
trajectory; while A^^cs varies wildly on a short time scale due to UV fluctuations, A^ is steady, 
and shows that the infrared fields never stray far from the vacuum. 

A^ will not meet all the conditions we set out for A'cs; for instance it will violate condition 
4. severely. However, it will be closer to a continuous function, since initial cooling removes 
UV excitation from the configuration, leaving weaker and more slowly varying fields for 
which the definition of E^Bf is less problematical. There is also less problem using coarsen- 
ing with this definition. For instance, on a 28^ grid, single coarsening after r = (5/4)a^ and 
double coarsening after r = 2.8(2a)^ and using an O(a^) improved definition of E^Bf, the 
discontinuity in A^ across the gradient flow separatrix is 0.9870 ± 0.0029 (drawing configura- 
tions from broken phase Yang-Mills Higgs theory at jSi = 7, see next section). To make A^ 
a continuous function modulo 1, the jump should have been 1. By rescaling A^ slightly, the 
discontinuity is removed almost altogether. The value of the discontinuity for larger lattices 
and deeper cooling is even closer to 1 with less noise. 

Let us check that the new definition of A^ will have A^ ~ 1/2 on the gradient flow 
separatrix. A point right on the separatrix will gradient flow to the saddlepoint configuration 
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Figure 9: Plot of J E ■ B from to r, as a function of r, for a series of points along a 
Hamiltonian trajectory as it passes through the gradient flow separatrix. 

and stick there.0 Perturbing the starting configuration slightly off the separatrix, it will cool 
to the saddle, miss slightly, and then slide off to a vacuum. The closer to the separatrix 
we start, the longer we will stick in the saddle before we slide out. There is some early r 
(transient) contribution to J J E'^Bfd^xdr while it is approaching the saddle, and then there 
is a contribution, almost exactly equal to 1/2, as it rolls from the saddle to the vacuum. 
By choosing tq large enough, we miss the transient and pick up only the 1/2, and so the 
N = 1/2 separatrix will correspond almost precisely with the gradient flow separatrix. We 
illustrate this with data from a Yang-Mills theory simulation in Fig. |^. The figure plots 



(shifted so the vacuum will have integer Nqs) against r for a series of points on a Hamiltonian 
trajectory as it goes through the separatrix. Each curve records the cooling process of a 
successive point on the Hamiltonian trajectory. We see that as the Hamiltonian trajectory 
approaches the (gradient) separatrix, the cooling path stays near the saddlepoint for longer 
and longer. When the Hamiltonian trajectory crosses the separatrix, the cooling paths roll 
out of the saddle towards the other side. The figure also shows the early r transient. We 
want to choose a tq large enough to eliminate this transient, so N ~ 1/2 will hold for a 
configuration which starts almost exactly on the gradient flow separatrix. 

^This saddlepoint is not the same as Klinkhamer and Manton's sphaleron, because we are considering the 
Yang-Mills Hamiltonian only, in a finite volume. But we know such a saddle will exist by the same argument 
Manton originally made for the existence of the sphaleron Gsl . 
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(22) 
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We should make sure that the result for Td will be independent of tq, once tq is large 
enough to kill the UV problems. We remind the reader that is computed by choosing an 
e -C 1, and computing the probability to be within e/2 of 1/2, 

r(l+e)/2 

P, = / PN(x)dx. (23) 

J{l-e)/2 

Also, one computes the time rate of change of A^, {\dN/dt\), evaluated for configurations 
with ~ 1/2. The rate is then 

rd = y^{\dN/dt\), (24) 

times the dynamical prefactor. 

Increasing the cooling time will decrease P^. The reason is that configurations which are 
near the saddlepoint at cooling time tq are spreading out from each other as cooling time 
progresses, as figure | illustrates. The rate of the spreading is given by the unstable frequency 
squared [uJ-)"^ of the saddle^. Increasing tq makes them spread further before we measure 
their Ncs', so the sample of states will be more diluted, by a factor of exp(— Aro/(ci;_)^). This 
reduces P^ by the same factor. However, we measure dN/dt by choosing two neighboring 
configurations on a Hamiltonian trajectory and finding the difference in their values of A^. 
The spread between these at cooling time tq will also increase as we increase tq, by the same 
amount; so {\dN/dt\) will go up by exp(Aro/(ciJ_)^). Hence, F^ will be Tq independent. This 
means, however, that neither {\dN/dt\) nor P/v(x) have simple physical interpretations. 

There is a modification of the above reasoning if tq is too short to eliminate the early 
transient; namely, a contribution to dN/dt due to the time evolution of the transient. In 
the complete calculation this will be compensated for because the dynamical prefactor will 
differ from the gradient flow value by a tq dependent amount, which becomes non- negligible 
at the same time the transient contribution to dN/dt does. In a complete calculation, F,;; 
will be independent of Tq, as we argued in the last section. 

We end this section by discussing briefly why we choose to define A^ based only on the 
Yang-Mills fields and using the Yang-Mills Hamiltonian for the cooling, rather than including 
the Higgs field. Doing so is reasonable because A'cs should be defined as a function of the 
gauge fields alone. Also, cooling all the fields under the full Hamiltonian is problematic, 
because the UV fluctuations of the gauge and Higgs fields renormalize the Higgs mass squared 
P6| . The bare potential needs a large, negative mass squared counterterm. However, the 
UV fluctuations are the first casualty of the cooling process, and so they stop generating 
a thermal Higgs mass squared, early in the cooling. To keep the minimum of the Higgs 
potential from changing radically, one would need to vary the bare Higgs potential in a 
complicated way as the cooling progressed. Other important fluctuation induced effects are 
also lost; for instance, the cubic term which makes the phase transition first order disappears 
as we cool the excitations. Depending on how we handle the Higgs potential during cooling, 
we will either cause symmetry to break during cooling when we start in the symmetric phase, 
or cause it to be restored when we start in the broken phase; we cannot avoid both, because 
there is a range of temperatures where each phase is metastable. Cooling only the gauge 
fields avoids this complication. 

^ which does not equal (tf_)^ of the Khnkhamer Manton sphaleron 
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4 Monte- Carlo calculations 



Here we discuss details of performing the calculation of Pn{x) and {\dN/dt\) using the 
definition of from Eq. ([T8|). We also discuss the computation of the dynamical prefactor, 
with and without "adding" hard thermal loops. Readers who don't do this kind of calculation 
will want to skip all but the first subsection and go to the results. 



4.1 multicanonical Monte-Carlo: the idea 

We need to compute the probability distribution of Chern-Simons number (really, the 
Chern-Simons number of a precooled configuration, defined in Eq. (|TB|)). The probability 
density that = x is 

PM = lim -J— [ VUVuV<^e-^^^^'''^^^Q{N{U) - x)e(x + Sx ~ N{U)) , (25) 

5x^0 ZOX J 

where U,u, $ are the SU(2) connection, the U(l) connection, and the Higgs field, and N{U) 
is defined in Eq. (|18|). Here / VUVuT>^ means the integral over the value of each field at 
each lattice site, and Z is the value of the integral without the step functions. On an 
lattice, this is a 16A^^ dimensional integral, which for = 40 is 1024000 dimensions. For 
this reason we turn to Monte-Carlo integration. In a canonical Monte-Carlo integration we 
generate a sample of configurations drawn with weight 

e-^^(^'"'*)pf/PuP$ , (26) 

and replace the integral with a sum over that sample. This will not do in the present context, 
because we want to know Pt<s{x) even where it is exponentially small. To get a good sampling 
there would require generating an exponentially large sample. 

We evade this problem by doing a multicanonical Monte-Carlo calculationH^. We rewrite 
Eq. (E5D as 



Sx^o ZSx J ^ 

(e-/(JV(c/))0(]V([/) - x)Q{x + 6x- N{U))) , (27) 

with f{x) some function we are free to choose. Now we generate a sample of configurations 
drawn with weight 

and replace the integration in Eq. (^) with a sum over this sample, with the term in the 
second set of parenthesis as the argument of the sum. By choosing /(x) ~ — lnPN(a;), we 
push the exponential suppression from the sampling into the integrand. The quality of the 
integration is now limited by how quickly we can generate a quality sample with this weight, 
and how well we can choose f{x), which we must do by some form of bootstrapping. 

A usual way to generate a canonical ensemble is by some Markov process. Given a 
configuration Ci and a realization ^ of some random noise distribution, the process returns 
a new configuration C2 = M(Ci,^). Define the probability to return a particular C2 as 

PiCuC2) ^ I di5{M{C,,i),C2). (29) 
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If M satisfies detailed balance, 



eMP{H{C,) - H{C2))) , (30) 



then iterating the Markov process generates the canonical distribution. We get the best 
statistics if applying M is numerically cheap and if the return value differs as much as 
possible from the starting value. 

If we have a Markov process which generates the canonical ensemble, we make it generate 
the multicanonical ensemble with weight function /(C) by making the following modification. 
Given a configuration Ci, generate C2 = M(Ci,^). Accept C2 as the next configuration in 
the sequence with probability 

min(l,exp(/(C2)-/(Ci))) (31) 

and otherwise reject it and make Ci the next configuration. This changes the detailed 
balance relation for the sequence to incorporate / into the weight. We may no longer want 
M(Ci, ^) to differ from Ci by as much as possible, though, because that may make the reject 
rate very large. Instead we want |/(M(Ci, ,^)) — /(Ci)| ~ 1. 



4.2 note on algorithm 

We need to do two kinds of things. The first is to evaluate the path integral for dimensionally 
reduced 3-D Yang-Mills Higgs theory. The second is to study the dynamics of the 3+1 
dimensional, classical theory. We comment briefly on the connection between the two; in 
particular we should compare the partition function of the classical theory to the path integral 



for the dimensionally reduced theory, a comparison first made by Ambj0rn and Krasnitz [|19 

The partition function of classical, 3-1-1 dimensional Yang-Mills Higgs theory (not worry- 
ing about the difference between lattice and continuum, which will not be important here) 
looks like 

Z = J VAiV^VEiVm{{D,E,y + (^/2)(inV''$ + h.c.)) exp{-H{A, $, E, U)/T) (32) 
H = jd'^{^ + ^ + (A$)^ + ) . (33) 

where E, the electric field, is the conjugate momentum of A, and 11 is the conjugate mo- 
mentum of $. The delta function enforces Gauss' Law. We can implement Gauss' Law by 
means of an adjoint valued Lagrange multiplier Aq, giving |T9| 

Z = j VA{D<^VEj:>IWAoe-KY>{iAl{{DiEiY + (^/2)nV'^<I) + h.c.)/T) exp(-//(A, $, E,I\)) . 

(34) 

The integrals over E and 11 are rendered Gaussian; doing them gives 

Z = JvAiV^VAoexp{-H{A,^,Ao)/T), (35) 

H = Jd^x(^ + + (A$)^ + . (36) 
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So the thermodynamics of the classical theory with Gauss' Law is governed by the 3-D 
dimensionally reduced path integral, but including the Aq field at zero (bare) Debye mass. 
We could get the path integral without the Aq field if we did not enforce Gauss' Law. There 
are two choices; either we treat the dynamics dropping Gauss' Law, or we include the Aq 
field in the 3-D Monte-Carlo parts of the calculation. 

The conservative approach is to include Gauss' Law in the dynamics, which means includ- 
ing an Aq field in the thermodynamics. The easiest canonical Monte-Carlo Markov process in 
this case is a short Hamiltonian trajectory starting with randomized but constraint respect- 
ing E and LI fields, with the multicanonical accept reject steps inserted between evolutions. 
This is a "constrained molecular dynamics" algorithm, and the problem of drawing E and 
n from the thermal distribution respecting the constraints is addressed in Alternately 



we could use heat bath and overrelaxation updates on the system described in Eq. (|35|) . 
Neither approach is terribly efficient. 

The other option is to assume that Gauss' Law is not very important to the dynamics, and 
not enforce it when we draw momenta E and LI from the thermal ensemble. This assumption 
is certainly justified with regards to measuring {\dN/dt\). What the combination of E and 
n forced zero by Gauss' Law would do if we did not set them zero is to gauge rotate the 
fields, but not the momenta. On long time scales that might be important, but at leading 
order in the length of a short Hamiltonian trajectory it only changes the gauge of the final 
configuration. Since is a gauge invariant object, this does not matter. A less rigorous but 
more physically intuitive way to see the unimportance of Gauss' Law to {\dN/dt\) is to note 
that it is roughly the instantaneous value of / d^xE^BI. Now the magnetic field is transverse 
by the Bianchi identity, so only the transverse components of the electric field contribute to 
dN/ dt. But Gauss' Law only depends on the longitudinal components, so dN/ dt should be 
the same whether or not we enforce it. 

It is less clear whether Gauss' Law will have a role in setting the dynamical prefactor, 
since it depends on longer time dynamics; but we are getting the dynamical prefactor wrong 
anyway if we don't enlarge the system somehow to account for hard thermal loop effects 
properly. We should deal with these two questions together. 

The chief advantage of not enforcing Gauss' Law is that there are very efficient update 
algorithms for the path integral without either the E, U or Aq fields, for instance the one 
developed by Rummukainen et. al. ||^. We adopt their lattice action, which is the standard 
Wilson 3-D Yang-Mills Higgs action, but we add a noncompact U(l) field. We use their 
update, which must be extended to include the noncompact U(l) field. In the noncompact 
formulation, the U(l) gauge field on a link is represented by a single real number Bi{x), and 
the terms in the action which depend on the U(l) field are 



-a 



^^(j)'^{x)Ui{x)exp{iaBi{x))(f){x + i) , (37) 



where a is the lattice spacing, z = tsua^Qw, and Ui{x) is the SU(2) connection on the i,x 
link. The compact formalism is the same except that {1/2){J2 B)"^ is replaced by a~^(l — 
cos{aJ2 B)). In the compact case the energy, as a function of one B, is a trigonometric 
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function. It is easy to perform an exact heat bath or overrelaxation step ||T0[- For the 
noncompact case the energy is the sum of a quadratic and a trigonometric function and it 
is not easy to perform an exact overrelaxation or heat bath update. We perform the update 
based just on the (much larger) quadratic term and include the trigonometric term by an 
accept reject step. The accept rate is quite high, so the cost to algorithm efficiency is low. 
We also occasionally gauge transform to bring the B fields towards Coulomb gauge so that 
the typical B is close to zero and the series expansion of exp{iaBi{x)) converges quickly. 

We always apply 0(a) improvement to the lattice action as described in |^ When- 
ever we refer to physical units in this paper they are always related to the lattice ones through 
0(a) improved relations. The improvement is essential to achieving the numerical accuracy 
we want at reasonable lattice spacings, and it makes small spacing extrapolations of most 
quantities unnecessary. For instance, we have computed the jump in 0^ at the equilibrium 
temperature for a = i/Tg'^T {(3 = 7) and a = {P = 10) to see how sensitive results 

are to varying a. The results are (A0^/(yf^T^) = 2.51 ± .02 and 2.56 ± .03 respectively; the 
(0(a^)) errors are smaller than the statistical errors we will be able to achieve for F^, so 
lattice spacing errors are under control. 



We have also used the technique of |2^ to study the infiuence of hard thermal loops on 
the dynamical prefactor. The idea is to add a large number of weakly interacting, ballistic 
charged particles to the lattice system, which reproduce the effects of the hard degrees 
of freedom left out when we set up our lattice. We refer the reader to for details. 
This approach demands that we apply Gauss' Law, and we must use the (Hamiltonian) 
thermalization algorithm as our Markov process. This is very inefficient, so we can only use 
the technique to study the dynamical prefactor, not the flux through the separatrix. 

4.3 finding Tc 

At what temperature should we study the sphaleron rate in the broken phase? We choose 
the equilibrium temperature for the phase transition. This is appropriate if the latent heat 
liberated during the cosmological electroweak phase transition is sufficient to reheat the 
universe to T^, which would be the case if the latent heat were large or the supercooling 
were small. It is in fact not clear whether this will be the case, in the MSM or in more 
viable extensions. We will choose Tc anyway, because we then have a well specified question 
with only one free parameter X/g"^; besides, complete reheating may generically occur in 
super symmetric extensions with a light stop. (It also may not; see [pO] .) 



In the effective 3-D theory we actually find the critical Higgs mass squared, which is 
related through the dimensional reduction procedure to the critical temperature. Since this 
involves comparing the thermodynamical favorability of the broken and symmetric phases, 
it necessarily involves some multicanonical technique. We will use a particularly simple 
approach, somewhat similar to the one used in |5^. It is based on the fact that, in a very 



long rectangular box, values of / ^"^^/V = {cfP') intermediate between the symmetric and 
broken phase values are obtained by having a mixed phase configuration, where part of the 
volume is in the symmetric phase and the rest is in the broken phase. 

The free energy, as a function of (0^), is — Tln(P), with P the probability to have that 
value of (0^). In the range of intermediate (0^), the free energy will vary linearly with (0^), 
since a change of (0^) represents a change of how much bulk free energy comes from one 
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phase and how much comes from the other phase. The slope of the hnear relation tells the 
free energy difference between the two bulk phases. This linear regime breaks down where 
(0^) comes close to the value in one or the other phase, since the phase interfaces then get 
close enough together to interact. 

Our approach is to add a (0^) dependent contribution to the Higgs mass, which we 
achieve by adding to the Hamiltonian a new nonlocal term. 



Choosing rj to be positive means that, if most of space is in the broken phase, the Higgs mass 
is heavier and the symmetric phase is favorable, whereas if most of space is in the symmetric 
phase, then the Higgs mass is smaller and the broken phase is favorable. The free energy is 
then a quadratic function, and the effective Higgs mass squared at its minimum, including 
the contribution from the rj term, gives the equilibrium Higgs mass. The added term is 
simple enough that we can modify the canonical update to incorporate it, and Monte-Carlo 
evolution will then naturally settle in a mixed phase configuration whose value of (0^) tells 
us the equilibrium Higgs mass. 



This approach can be viewed as a type of multicanonical Monte-Carlo with the rj term as 
the multicanonical reweighting. 

In practice we start with a long but very narrow box and a high value of rj, to get a 
preliminary value. The narrowness is necessary to make it easy to nucleate a bubble of one 
phase in the other. To get the large volume limit, though, we need to go to a wider box; we 
necessarily need results in a regime where one phase cannot easily nucleate in the other. We 
get an initial condition for a box an integer number of times wider in each short direction 
by extending the final configuration in the skinny box periodically. We also use a smaller 
value of T) to improve the resolution of the determined m^. Our final values for are 
typically taken with a box of dimension, in physical units, of x {96 /g'^T), easily 

large enough to achieve the large volume regime. 

4.4 practicing with the symmetric phase in small volume 

Before presenting the determination of the broken phase sphaleron rate, we will do a "practice 
run" on a system where we can get good statistics more quickly, which is the symmetric phase 
in a small enough volume to suppress topology change. This problem is almost certainly of 
no phenomenological significance, but it will let us test the tq dependence of our technique 
and to study whether hard thermal loops do indeed modify the dynamical prefactor. 

For the time being we drop the U(l) subgroup. We choose a very weak scalar self-coupling 
oi X = 0.025, and a large Higgs mass squared (m^ ^ l.bg'^T'^), so we will be firmly in the 
symmetric phase. We use a 12^ lattice with a lattice spacing of a = l/{4g'^T) {j3g = 16), 
so the physical volume is {3/g'^T)^. This volume is small enough that topology change is 
suppressed and the broken phase techniques we are developing are applicable, but not so 
small that it will be hard to gather good statistics. It is also small enough that numerical 




(38) 




(39) 
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N using Tq=1.25, 6.25 

Figure 10: Probability distribution of for two values of tq: tq = 1.2'ba? (the rounder 
curve which is lower at = 1/2); and tq = 6.25a^ (the curve which is higher at = 1/2). 

costs are not overburdening, so we will enforce Gauss' Law (and use the less efficient update 
algorithm) . 

We measure N by cooling for r = (5/4)a^, coarsening once, and using an 0{o?) improved 
definition of EfBf during the subsequent cooling. The first thing we do is to determine the 
actual discontinuity in across the gradient fiow separatrix, which will differ from 1 because 
of lattice artifacts. To do this we generate an ensemble of configurations on the gradient 
fiow separatrix. For each, we perform a short Hamiltonian evolution which crosses the 
separatrix, and we measure A^ at closely spaced intervals during the crossing to determine 
the discontinuity; it is AA^ = 0.982 ± 0.005, where the error is the standard deviation[^ 
This tells us how to rescale all further measurements of A^ so that there will be an integer 
discontinuity at the separatrix (up to acceptably small noise). 

Next, we measure the probability distribution for A^. We only need to do this in the range 
0<A^<l/2, by periodicity; so we add 1 to all negative values of A^ and then put N > 1/2 
into range by A^ — 1 — A^. We use a continuous, piecewise linear reweighting function, with 
the widths of the linear pieces chosen by hand and the slopes of each segment determined 
by an automated bootstrapping procedure. We computed the probability distribution for 
two values of tq; tq = 1.25a^ and tq = 6.25a^. For the tq = 1.25 data we also recorded 
N{to = 6.25); the data can therefore be used to get the probability distribution for either 

^°It is important to compare the extrapolation based on a fit of the last few points on one side of the 
discontinuity with the value on the other, to remove errors from the time step size of the Hamiltonian 
evolution used to find the discontinuity. 
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Figure 11: A winding number changing section of a Hamiltonian trajectory in constrained 
volume Yang-Mills theory. After the crossing the system settles immediately into the neigh- 
borhood of a topological vacuum. 



choice of tq. 



The probability distributions are compared in Figure |T0 
different 



The distributions are clearly 
Cooling longer before integrating E'^B'^ concentrates probability around = 



and thins out the large A^ configurations. It is less likely to have A^ within some range of 1/2 
for the larger tq. But this does not mean that the different tq values give different values for 
the diffusion constant, as we still have not included {\dN/dt\) or the dynamical prefactor. 

To get the dynamical prefactor and {\dN/dt\) we first need a sample of configurations 
very close to A^ = 1/2. We get them by multicanonically sampling, not necessarily with 
the same reweighting function used to find the probability distribution of A^. Then, for 
each, we choose momenta out of the appropriate distribution and perform a Hamiltonian 
evolution, with the algorithm of [|18]. Once the Hamiltonian trajectory has settled into the 
neighborhood of a vacuum, we return to the configuration before we started the Hamiltonian 
evolution and reverse the sign of the momentum; then we evolve. Since momenta are odd 
and fields are even under time reversal, this computes the Hamiltonian trajectory in the 
backwards time direction. We determine {\dN/dt\) from the first time step in each time 
direction, and the dynamical prefactor from the number of A^ = 1/2 crossings, as discussed 
in the previous section. It is also easy to find how many times the Hamiltonian trajectory 
crosses the gradient flow separatrix, since / E°'B^d^xdT abruptly changes sign when it does; 
so we can also identify how well correlated crossings of the A^ = 1/2 and gradient flow 
separatrices are. (We cannot directly determine the gradient flow separatrix prefactor if the 
correlation is not good, because our sampling procedure is for the N = 1/2 separatrix and 
not the gradient flow one.) 



Figure O gives an example of a Hamiltonian trajectory developed in this way, with 



To = 3.75a^. It is clear in this figure that once the trajectory settles in the neighborhood of 
a topological vacuum, in the sense that A^ ~ 0, then it stays there for some time. Since the 
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v^ilnp Pit Tn = 1 


vpilnp p^t Tn = 

V CAil HVj CA U I I 1 • ^ t_/ L(j 


InPflA^ — 51 < 05) 


—13 39 ± 22 


— 14 30 ± 18 


\dNldt\ 


{ 121 ± 004')o^T 


{ 174 ± 005)o^T 


prefactor 


.52± .03 


.554 ± .024 


ln(r/(a^T)4) 


-7.03 ± .23 


-7.51 ± .19 


r 


(.0009 ± .0002)a^T4 


(.0005± .0001)atT4 



Table 2: Ingredients and results for the A^^cs diffusion rate in a cubic, periodic volume 'ij g^T 
on a side. The measurements with two values of tq disagree at 1.6a. 



maximal Lyapunov exponent of the Yang-Mills Higgs system is known to be about 0.3(7^T 
|5^ , |19|, the direction of the next permanent change in A^cs will surely be statistically 
independent from the previous one. This fact is essential to the whole approach; our most 
basic assumption is that the very long time A^cs diffusion is made up of a series of statistically 
independent steps, and we only need find out how frequently one of those steps is taken 
(which is what we are computing when we say we are computing 7^). 

Also note from Figure |ll| that there is no "overshoot" after reaching A^ = 0, no sign 
that the trajectory is continuing in the direction of the next separatrix. This is true of all 
trajectories we have studied, both in finite volume and in the broken electroweak phase, 
which confirms the absence of prompt multiple crossings. 

Our final results for this small volume system are presented in Table |^. Note that both 
the probability to be near the separatrix and (|(iA^/(it|) vary quite a bit when we change tq, 
but in opposite directions. Also, the shorter cooling leads to a smaller dynamical prefactor, 
though for the volume and cooling considered here the difference does not turn out to be 
large. The determined value of F for the two values of tq differ at l.Gcx. This does not bother 
the author since it is the first statistical fluctuation above 1.5ct he has encountered since 
starting numerical work; he had one coming. 

We can make a more stringent check of the tq independence of the method by using the 
data set taken with tq = 1.25a^ to determine the probability distribution of A^ for each value 
of To, since we recorded N(tq = 6.25a^) at each point while developing this data set. This 
data set gives a slightly different probability to be at large N{t = 6.25a^); the —14.30 in 
the table becomes —13.88. Using this number, we get ln(F/(Q;^T)^) = —7.09. To compare 
to the To = 1.25a^ data we must remember that the probability distributions are now 100% 
correlated, so the expected difference is just from statistical errors in {\dN/dt\) and the 
prefactor. The difference in the logs of the rates is .06 ± .09, so they do agree within error. 
The results for Yd are indeed tq independent. 

How do hard thermal loops change the rate? At the level of thermodynamics, hard 
thermal loops become just a Debye mass. They make the Aq field heavier, which pushes it 
further to the regime where it decouples. So their thermodynamic influence is very small. 
Hence, they will barely change the probability distribution of A^, which will have a well 
behaved large HTL strength hmit. Similarly, as we have discussed, {\dN/dt\) does not 
depend on longitudinal E fields, and hence depends on the Debye mass only through its 
thermodynamic influence on the gauge fields. The flux through the separatrix should depend 
weakly on hard thermal loops and should have a good large HTL limit. 
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20 40 60 80 

time xg^T 

Figure 12: Hamiltonian trajectory through the separatrix in the symmetric phase in small 
volume, including hard thermal loops with m|, ~ ASg^T"^. Plasma oscillations are clear, and 
there are numerous crossings of the separatrix. 



The dynamical prefactor is a completely different matter. It is a dynamical quantity 
which depends on unequal time correlators, potentially over quite long times. Hard thermal 
loops will change the time evolution of infrared degrees of freedom on all time scales longer 
than the inverse plasma frequency. Arnold, Son, and Yaffe (ASY) have argued that hard 
thermal loops will suppress the baryon number violation rate by a factor parametrically of 
order (g^T^/mj^), because the number of crossings of the separatrix per permanent A^^cs 
change will be of order mj^/g^T'^. That is, they predict that turning on hard thermal loops 
will reduce the dynamical prefactor by an amount linear in mjy when it is large. Their 



arguments have recently been verified in the symmetric phase [^. Now we need to check 
what the hard thermal loops do in the broken phase. 

To address this question we include the hard thermal loops by using the technique de- 



veloped in |26[. We add a large number of new "particle" degrees of freedom to the system, 
which propagate the hard thermal loop effects. For now we will be happy to know what the 
hard thermal loops do when they are extremely strong, so we can understand the parametric 
limit in which the ASY arguments should hold. With this in mind we put in particles of 
charge Q = 0.1 and number density 50/a^,[3 which give a Debye mass of 43(7^T^, an enor- 
mous number about 10 times the value in the MSM at a realistic value of g^. A Hamiltonian 



trajectory crossing the separatrix is shown in Fig. 12. The qualitative features are indeed 



the same as Arnold, Son, and Yaffe predict, see Fig. 8 of plasma oscillations drive the 

^^See for the implementation and the definitions of these quantities, and their relation to the Debye 
mass. 
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system across the separatrix numerous times. The dynamical prefactor is correspondingly 
significantly smaller than without hard thermal loops; for these parameters it is about 0.16. 
Note however that we have had to add truly huge HTL effects to achieve this value, so at 
the realistic value the effect may not be too significant. We will study this question for the 
broken phase case in the next section. 

How do our nonperturbative results compare to perturbation theory? We will not attempt 
to do a complete one loop calculation of the sphaleron rate in finite volume, but it is quite 
easy to compute the "sphaleron" energy, the energy of the saddle point between topological 



minima. We can use the technique of |£4[, or any other technique which can find a saddle 
point solution (we have one). We find E = 27.77/A^, with the linear dimension of the 
lattice. In our case that means (3E = 37 and we would naively expect a rate suppression of 
exp(— 37) before zero modes and the fiuctuation determinant are included. This might be 
compared with the free energy difference between = and N = 1/2, which is of order 
15. We know that the complete inclusion of the zero modes and the fluctuation determinant 
is likely to make up some of this difference, but it certainly will not account for all of it. 
Nonperturbative physics is at work and it enhances the rate of sphaleron transitions in this 
case. 



4.5 broken phase rate 

Now we will apply the same technology to the physically interesting case of large volume, 
broken phase SU(2)xU(l) Higgs theory. Since the previous subsections already explained 
both the Monte-Carlo update technique and the real time tools used to flnd {\dN/dt\) and 
the dynamical prefactor, we will just discuss here how this case is different from the small 
volume, symmetric phase calculation, and what we have to do differently to get it to work. 

As we discussed, we will work at the critical temperature, which corresponds in the 3-D 
language to the critical Higgs mass. This means that the broken phase, which we want 
to study, is actually only metastable; the symmetric phase is equally thermodynamically 
preferable. In a small volume there can be fairly easy tunneling between the two, but the 
metastability becomes stronger as the volume becomes larger. We must choose a volume 
which is large enough that metastability is very strong and tunneling between the phases 
will not occur. This means that the physical size of the lattice we use must be signiflcantly 
larger than the physical size of the sphaleron, which would have to be true anyway to keep 
the exponential tails of the sphaleron from seeing each other around the periodic boundary 
conditions. 

The problem of tunneling to the symmetric phase is made worse because the sphaleron 
has a zero of the Higgs condensate at its core, so it looks something like a symmetric phase 
bubble. To keep from nucleating to the symmetric phase we should use a volume big enough 
that the free energy of the state intermediate between phases is comparable to the free 
energy of the sphaleron. We have studied three values for x = X/g"^, x = 0.047, 0.039, and 
0.033; for the former two we used a physical volume of {16/ g'^TY and for the latter we used 
(13.33/5'^T)^. These were all sufficient to prevent nucleation of the symmetric phase, but for 
X = 0.047, a volume of {12.8/ g'^Ty was not. The volume requirement becomes less severe 
as the phase transition becomes stronger at smaller x, so for x = 0.039 and 0.033 we used a 
good margin of excess volume. The volume requirement would also have been less severe if 
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we worked below, rather than at, the equihbrium Higgs mass parameter. 

The need for a large volume drives up numerical costs in two ways. One is obvious; we 
need to update a lot of lattice volume which is "dead weight" since the sphaleron is not 
sitting there. But the large volume also makes the multicanonical algorithm perform worse. 
Examining figure 0, we see the free energy rises roughly linearly with at small values of 
A^. This behavior is also expected analytically in the broken phase case, see for instance [^|. 
Naively, then, having the gauge fields being part way up the sphaleron at one place in the 
box is not thermodynamically favored over having them half as far up the sphaleron in two 
different places, at least for relatively small N. In fact the entropy from getting to choose 
the locations of the two places means that this may be slightly preferred to being further up 
the sphaleron in one place, for the range of where the free energy is varying linearly with 
A^. But at larger A^ the free energy for a single "near sphaleron" levels off while that for two 
continues to rise linearly, so near N = 1/2 we prefer having a single sphaleron. Somewhere 
in between there is a mismatch in what kind of configuration is dominating the ensemble, 
and such a mismatch can reduce the efficiency of the Monte-Carlo. 

To cut the numerical demand we integrate out Aq fields (i.e., we do not enforce Gauss' 
Law when we study dynamics), which allows us to use the very efficient update algorithm 
of Rummukainen et. al. [Q. In fact a single step of the Rummukainen et. al. update 
(one heat bath and four overrelaxation sweeps) is far too large an update of the fields; if 
the value of f{N) typically changes by more than 1 under one update then the accept rate 
for updates becomes very small. Instead, we alternate between performing a "scaled back" 
version of the update in which only some "mod p checkerboard" of the sites are updated, 
and an overrelaxation sweep updating only the Higgs fields. Since the overrelaxation sweep 
of the Higgs fields does not change the gauge fields, and since A^ depends on the gauge fields 
alone, this sweep can be automatically accepted. Doing this means that the evolution of 
the Higgs fields can be made "fast" in comparison to the update of the gauge fields. That 
ensures that the thermalization of the gauge fields is not gummed up by slow evolution of the 
Higgs fields. Even with the extra updates for the Higgs fields, and the blocking procedure for 
accelerating the measurement of A^, most machine time is spent measuring A^, and further 
improvements in the update algorithm won't help. 

The results for the fiux through the N = 1/2 separatrix are given in Table |^, and also 
in [jl6|. We have also measured the dynamical prefactor for the x = 0.039 case, without 
added hard thermal loops or enforcement of Gauss' Law. The value we get is .33 ± .05, 
slightly lower than the value for the gradient fiow separatrix, which is 0.40 ± .05 (we measure 
this by choosing a value of Tq large enough that the two separatrices almost coincide, so 
there is a one to one correspondence between crossing A^ = 1/2 and crossing the gradient 
fiow separatrix). It is not clear whether this represents some interesting dynamical behavior 
of the theory or whether it means that the Yang-Mills gradient fiow separatrix is not the 
optimal divider between topological vacua. 

Another interesting question is how the dynamical prefactor depends on hard thermal 



loops. Unfortunately, the numerical cost of using the "particles" technique of |]26| is so high 
in this context that we have to cut a few more corners to make the calculation. We drop 
the U(l) factor, increase the lattice spacing from a = 2/{bg^T) to a = l/{2g'^T) {P = 8), 
and reduce the volume to {14/g'^T)^. To prevent nucleation to the symmetric phase, we 
work at a somewhat larger value of x, x = .042, and below the equilibrium temperature. 
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so the broken phase is more stable than the symmetric phase. We choose the temperature 
(the thermal Higgs mass, really) so the Higgs condensate is 0o = l-^S'T, just more than 
enough to make the sphaleron rate too low to erase baryon number (we expect). We cannot 
directly measure the sphaleron rate with this set of parameters because the update including 
particles is expensive enough that the multicanonical calculations are prohibitive, but we can 
get a good sample of points on the separatrix by using a reweighting which favors ~ 1/2 
strongly, and we can get sufficient statistics for the dynamical prefactor to make a good 
determination. 

We can "shut off" the hard thermal loop contribution to the dynamics without changing 
their contribution to thermodynamics by not allowing the particles to move during the 
Hamiltonian evolutions used to determine the dynamical prefactor. Also, if the arguments 
of Arnold are correct |5^, changing the velocity at which the particles move changes their 



contribution to the key part of the hard thermal loops, the ci; <C ~ g^T regime, linearly 
in the velocity. Hence we could explore very strong hard thermal loops either by putting in 
very many particles, or by making them move very fast. The numerical cost is the same but 
the memory costs favor the latter. 

We use a realistic total hard thermal loop strength, m|, = bg^T"^ from particles. (The 
standard model value is m|, = (11/6)(7^T^, and since g"^ ~ .40 on performing the dimensional 
reduction calculation [^, this is just smaller than the value we used.) The dynamical 
prefactor changes only mildly, from .52 ± .05 to .40 ± .05, when we turn the particles on. 
When we increase the particle velocities to 4 times the speed of light, the prefactor becomes 
.15 ± .03. Hard thermal loops do indeed reduce the dynamical prefactor, which is already 
less than 1 without them. However, the parametric limit in which the reduction is large is 
not achieved for realistic parameter values. Also note that the value .52 is larger than we 
got at X = .039 and = We assume this is because the larger value was evaluated 

below the equilibrium temperature, where the Higgs condensate is larger and stiffer and 
the sphaleron is smaller and more energetic; its decay should be more vigorous and less 
susceptible to buffeting by large IR fields. 



5 Results 

We present our results for the diffusion constant for Nqs in Table ||. The first quoted value is 
without the dynamical prefactor. The dynamical prefactor is less than 1 for the (A^ =1/2) 
separatrix we have used, and also for the gradient flow separatrix; its value for the x = 0.039 
data is about 0.33, and we will take this to be representative of the other two values of x as 
well. It is not clear whether the prefactor is less than 1 because our separatrix is sub-optimal, 
or because the dynamics are nontrivial in a way which often leads to multiple crossings. 

We have found that the dynamical prefactor depends on the strength of hard thermal 
loops; but for the realistic Debye mass the effect is not strong. A reasonable estimate for 
the relevant dynamical prefactor for the cases of interest and for the value of tq we used to 
deflne the measurable A^ is around 0.3 ± 0.1. We have included a row in the table where we 
include this (estimated) effect in the rate. 

For amusement, we have also compared the broken and symmetric phase probability 



distributions for A^. We present the results in Figure 13, which also shows how varies 
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Figure 13: Free energy (left) and (^Vbroken — 0Vsymm)/(5'^^^) as functions of at x = 
Xjg^ = 0.039, in a {16/ g'^T)^ volume (at tq = 3.6/{g'^T)^). In each case the upper curve is 
the broken phase and the lower curve is the symmetric phase. The plot of 0^ shows that 
the volume used was large enough that the sphaleron did not bring us anywhere near a 
transition to the symmetric phase. The behavior of the two phases is completely different; 
in the broken phase there is a free energy barrier, and in the symmetric phase there is not. 

with in each case. The two data sets were taken using identical values for all parameters 
(lattice spacing a = 2/5g^T, x = 0.039, critical temperature (i.e. critical Higgs mass), lattice 
volume = (IQ/g'^T)^, and tq = 3.6/ {g'^Ty), but starting with a broken phase initial condition 
in one case and a symmetric phase initial condition in the other. A barrier to changing A'cs 
is clearly present in the broken phase case and clearly absent in the symmetric phase case. 

We cannot use the data for the symmetric phase case to determine in the symmetric 
phase. Although it would be straightforward to compute {\dN/dt\) and get the flux through 
the separatrix, the calculation of the dynamical prefactor is impossible. A Hamiltonian 
path through the separatrix does not settle into the neighborhood of a topological vacuum 
and stay there for a long time; it just continues to wander around, as we already saw in 
Figure ^. But we can measure F^ in the symmetric phase with purely real time techniques; 
the value including hard thermal loops (and actually for pure SU(2) Yang-Mills theory) is 
Frf = (29 ± 6)al^T^, or — InF^ ^ 13.9. (This result may get revised downwards somewhat 
when the issues involving logarithmic corrections |^ have been fully accounted for.) As 
expected, the broken phase rate at small x = X/g"^ is enormously smaller; for x = .033 the 
ratio is about 10^. 

5.1 comparison to perturbation theory 

We want to compare the determined sphaleron rate to two things; perturbation theory, and 
the value required to avoid erasure of baryon number generated at the electroweak phase 
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transition. One loop perturbation theory gives 00 



7 



r, = 4T'^^(^^) j A/-.^v„,.e-^'=- . (40) 

Here 0o is the broken phase Higgs condensate expectation value, uj^ is the unstable frequency 
of the sphaleron, A/trA^Vrot are zero mode factors, k is the one loop fluctuation determinant, 
and Esph is the energy of the Klinkhamer Manton sphaleron, using the tree level Hamiltonian. 
For small X/g"^ = x, —Tlnn equals the energy due to the one loop effective potential term, 
plus a modest correction [|1^]. We can guess that the dominant two loop corrections to Eq. 
(^) are absorbed by including the two loop effective potential terms in the Hamiltonian. So 
it seems reasonable to estimate the sphaleron rate by Eq. (|^), but setting k, = 1 and solving 
for the sphaleron energy using the two loop effective potential at the equilibrium temperature. 
One should also solve for the zero modes and uj-/(j)o at this value, but they are very weak 
functions of the effective potential [0. We use the values from at a; = (X/g^) = 0.04 



for these, but solve for the sphaleron energy, Esph = 47ri?0o/5'! numerically, using the two 
loop effective potential at Tc. We use the two loop potential presented in |5^, without pieces 



from longitudinal gauge bosons (assumed integrated out). We also drop two loop terms 
proportional to Xg^ or A^, because the perturbative determination of 0o is an expansion 
in X/g"^, and such terms contribute at the same or higher order as unknown 3 loop terms. 
(Including these 2 loop terms moves 0o closer to the nonperturbative value by an insignificant 
amount.) The "two loop" analytic sphaleron rate, also included in Table 0, is about exp(2.5) 
times faster than the numerically determined nonperturbative rate, and falls further off when 
we include the dynamical prefactor^ The difference is more than can be explained by the 
difference in (pQ, but it is not huge in the sense that it represents a change of less than 20% in 
the exponent. The difference gets smaller, relative to the exponent, as the sphaleron energy 
becomes larger. 



5.2 comparison to the erasure bound 

We should compare the sphaleron rate to the limit set by requiring that baryon number not 
be erased. The rate at which sphalerons degrade baryon number is 

1 rfiVB _ 13iVp^ . 



Nb dt 



^dT-^ (41) 



where iVp = 3 is the number of families, and the numerical factor 13A'^p/4 would be smaller 
in theories, such as supersymmetry, in which additional degrees of freedom can store baryon 
number.[^ Integrating from the end of the phase transition to the present day, 

ln(iVB/iVB(T,)) = / T,{T{t))T-\t)dt, (42) 

4 Jta 

^^The definition of T used in ^ is the response rate to a chemical potential, which is half the diffusion 
rate |3^; so Eq. (|4C| ) differs by a factor of 2 from the expressions in those references. 

^•^There is no literature calculation of the dynamical prefactor, including hard thermal loops but in the 
perturbative context. 

Again, there is a factor of 2 difference from the reference because they write in terms of the response to 
a chemical potential, which is half the diffusion constant. This 2 cancels the other 2. 
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where we have shown the dependence of Ta on T and of T on t. 



Now InTii is very sensitive to y {<f)^(f)), and hence to T; so we can approximate Inrrf(T) 
InTdiTc) + (T — Tc){d\nTd/dT)\T=Tc, and perform the integrah 



ln(-ln(iV3/iVB(Tj)) =ln(39/4) + lnf^#l') -In^ rflnr.(T(t)) 



T4 / \ Trft 



(43) 



By the chain rule, 

dlnVd dlnTd dy dliaT 
Tdt ~ dy dT dt ' ^ ^ 

where y = m'^{T)/{g^T'^) is the dimensionless thermal Higgs mass squared. 
We get dy/dT from the 1 loop correction to ||57|| , 



dy^ 8X + 4g^^ + g\3 + tan^ew) , . 

dT~ Sg^T ' ^ ^ 

and we get dlnT/dt from the Friedmann equation in a radiation dominated universe, 



At"^ 3 30 dt ^ 45 mpi ' ^ ^ 

where g^ is the number of radiative degrees of freedom in the universe {g^, = 106.75 in 
the minimal standard model) and rupi ~ 1.22 x lO^^GeV is the Planck mass. Finally, we 
determine dlnTa/dy perturbatively, by varying y slightly from the equilibrium value and 
recomputing the two loop sphaleron rate. The dependence is quite strong. We include it in 
Table 

The most widely cited discussion of baryon number erasure after the phase transition 
makes the approximation that the baryon number violation rate after the phase transition 
is constant for about one Hubble time |58|. In fact, because depends very strongly on y, 



which in turn depends strongly on T, most baryon number erasure occurs in the first 10"'^ 
Hubble times after the phase transition. Hence the initial rate of baryon number violation, 
Frf(Tc), which prevents washout is 10^ times larger than assumed in |5^, leading to a weaker 
bound on Fd(Tc), roughly 

- ln(Frf(Tc)T-^) > 30.4 - ln(Te/100GeV) . (47) 

The values of g^, and dy/dT will both be larger in supersymmetric extensions of the standard 
model. Unless quite a number of supersymmetric partners have masses under 100 GeV, which 
now seems unlikely, g^, will not be too much larger; but if there are stops with SUSY breaking 
masses of less than 100 GeV, which is necessary to get a strong enough transition without 
violating the current experimental Higgs mass bound, then T dy/dT gets extra contributions 
from stops which bring it up by about a factor of 2. The bound, Eq. (|47|), is weakened by 
about 1. Also note that, because Eq. (^) is for the double log of Nb/Nb{Tc), failing to meet 
the bound by 1 means the baryon number is diminished by exp(exp(l)) ~ 15, and failure by 
2 reduces baryon number by exp(exp(2)) ~ 1600; so the bound is quite sharp. 
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6 Conclusion 



We have shown how to define a lattice measurable which allows a (multicanonical Monte- 
Carlo) nonperturbative determination of the broken phase sphaleron barrier height. We have 
combined this with real time techniques to measure the diffusion constant of A^^cs (and hence 
baryon number violation) in the broken electroweak phase nonperturbatively, including the 
dynamical prefactor. We find that the diffusion constant is smaller by about exp(— 3.6) 
than in a perturbative estimate using the two loop effective potential and no wave function 
corrections (and assuming a dynamical prefactor of 1). The difference is too large to ascribe 
to the difference between the perturbative and nonperturbative value of 0o- However, it 
represents a shift in the exponent of less than 20% from the perturbative estimate. 

We have also demonstrated that the physics of hard thermal loops does change the 
sphaleron rate in the broken phase, apparently consistently with the arguments of Arnold, 
Son, and Yaffe |^3|]. But to really achieve the parametric limit they discuss takes an unreal- 
istically large Debye mass; for physical values of the parameters, the correction due to hard 
thermal loops is fairly minor. 

Interpolating between the values oi x = \/ where we have measured, and including 
the estimate for the dynamical prefactor, we get a bound of about (X/g^) = x = .037 in 
the standard model and x = .039 in the MSSM when it can be perturbatively reduced to a 
standard model like effective theory. These are slightly looser than we quote in |]16[ because 
the measured value of the dynamical prefactor is smaller than our estimate there. Another 
convenient way to state our result is that the bound on the Higgs condensate after the phase 
transition is about 0o = l.G7gT in the MSM and 0o = l.eO^fT in the MSSM (with the MSSM 
value weaker because the larger temperature dependence of the thermal Higgs mass makes 
the erasure rate fall off faster with time after the transition). 

The bound is softened if the universe does not reheat to Tc during the phase transition, 
because the phase transition then ends at a lower temperature with a larger Higgs vev. 
Incomplete reheating may well be generic. Unfortunately, we do not have a nonperturbative 
measurement of the bubble nucleation action, which is required to determine definitively 
whether reheating happens. This is an interesting project which can perhaps be approached 
by techniques similar to what we have used here, i.e. defining a separatrix corresponding 
to the critical bubble (as the separatrix here corresponds to the sphaleron), and making a 
complete calculation including the dynamical prefactor. 

We should comment on what we expect to be true below the equilibrium temperature 
and in the case of a light stop. As the temperature falls below equilibrium and the Higgs 
condensate becomes larger, perturbation theory should become better, at least in the sense 
that the error in the exponent should get smaller compared to the magnitude of the exponent. 
It would be very surprising if the nonperturbative rate switches to being faster than the 
perturbative estimate, though. It is straightforward, though expensive, to repeat the analysis 
here for temperatures below the equilibrium point. Probably it will only be worth it after 
we have nonperturbative information on the bubble nucleation rate. In the case of the 
light stop, we expect the most important differences to be in the effective potential, in 
which case perturbation theory should do as well as it does here, if we permit a "by hand" 
correction for the strength of the phase transition. We expect this because, while the stop 
contributes at one loop to the effective potential, scalar interactions only influence wave 
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function corrections (say, the Higgs gradient energy of the sphaleron) at two loops, and 
there is no direct interaction between a hght right stop and the SU(2) gauge fields. 
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